Non Classical Rotational Inertia Fraction in a One Dimensional Model of Supersolid 
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We study the rotational inertia of a model of supersolid in the frame of the mean field Gross- 
Pitaevskii theory in one space dimension. We discuss the ground state of the model and the existence 
of a non classical inertia (NCRI) under rotation that models an annular geometry. An explicit 
formula for the NCRI is deduced. It depends on the density profil of the ground state, in full 
agreement with former theories. We compare the NCRI computed through this theory with direct 
numerical simulations of rotating ID systems. 



INTRODUCTION 

Since pioneering works by Andreev and Lifshitz[T], 
Chester [5], Leggett[3] and others, supersolids have been 
dreamt up as a kind of Bose-Einstein condensation of 
defects, vacancies or interstitials. They would achieve a 
coherent state that could allow a matter flow trough the 
crystal. Although the quest for a supersolid state over the 
last 40 years has failed[4|, the context has totally changed 
with the recent experiments by Kim and Chan|3{7]. In 
these experiments, solid Helium^ fills a torsional oscilla- 
tor under small oscillations and the rotational frequency 
is measured. Surprisingly, the rotational inertia shows a 
drop below few tenths of Kelvin. This non-classical ro- 
tational inertia (NCRI) is believed to be the signature 
of the transition of a fraction of the solid into a super- 
solid state. The situation has become puzzling as other 
experiments have been performed. Thus, although the 
drop of the moment of inertia has been confirmed, crys- 
tal annealing was shown to lower dramatically the am- 
plitude of this NCRI 8, 9J. Similarly, solid Helium sub- 
mitted to pressure gradient could not flow except when 
large grain boundaries where present in the sample |10l - 
[T^ . Moreover, the responses of solid ^He to a localized 
pressure jump presented no evidence of superflow in the 
solid[T31[Il] .The experimental context presents thus ap- 
parent contradictions between NCRI measurements and 
pressure driven flows with the role of disorder (through 
vacancies, grain boundaries...) to be elucidated. On the 
other hand the theoretical framework for describing su- 
persolids presents also fundamental puzzles (see the re- 
cent review of Prokof 'ev[TS] where the influence of the 
disorder is particularly discussed). Beside Penrose and 
Onsager argument [16] , Monte-Carlo models claim that a 
perfect crystal cannot exhibit supersolid behavior pTllTB] . 
However, the account for exchange processes between 
neighboring atoms [3 [TO], the densities and the role of the 
vacancies in the dynamics raise still fundamental ques- 
tions on the existence and the nature of the supersolidity 
(see [101 m] for instance). 

An alternative issue consists of using the Gross- 



Pitaevskii (CP) model[55j to describe the dynamics of 
a quantum solid, as proposed in 1994 by Pomeau and 
Rica|23j. The original GP equation [22] is used, with 
a roton minimun in the dispersion relation, where the 
ground state exhibits a first order phase transition to a 
crystalline state. However, the assumptions underlying 
the GP equation are not strictly speaking valid for He- 
lium although this equation is believed to give a good 
qualitative description of superfluid Helium. Therefore 
this model, even crude in its basic structure, is at least 
a good testbed for theories of supersolids, that are still 
in a state of uncertainty. In Ref. [M] US] , two of us have 
developped the theory for the long wave perturbations 
of this model of supersolid and we have shown that this 
model was able to conciliate the apparent experimental 
contradiction discussed above. In the present paper, we 
study the one-dimensional (ID) version of this model. 
Beside the simplicity of the ID approach, which then al- 
lows precise determination of the crystal structure, the 
ID limit is particularly interesting since it can model to 
some extent an annular geometry. 

THE MODEL 

The starting point is the original GP equation[22] for 
the complex wavefunction ^{x,t) in one space dimension: 

where U{s) is the two body potential depending on the 
relative distance. The potential U{s) should satisfy 

/oo 
U{s)ds < 00, 
-oo 

for stability and its Fourier transform: 

/oo 
U{s)e''''ds (2) 
-00 

has to be bounded for all k. Moreover, as we will see 
later, we shall require also that the Fourier transform Uk 
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becomes negative at some kc to allow roton crystalliza- 
tion. 

This dynamics conserves the Hamiltonian (or the en- 
ergy, following dtij ^ 4§,) 



H 



2m 
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the right values of the speed of sound c, and the three 
roton parameters 1271. 
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With x' — X / a, t' — and i/)' ~ y'noV'i the dimen- 
sionless GP equation for the Heaviside interaction (we 
drop the primes hereafter) reads: 



Ui\x-y\my)fmx)\'dydx, 



CO ^ — OO 



the number of particles 
N = 



\'ijj{x)\^dx 



and the linear momentum 
ih f°° 

P = --/ ird,^l>-i;d^4,*)dx. 

^ J -oa 

According to the energy, the ground-state solution is real 
since any non uniform phase increases its energy. 

The dynamics exhibits indeed an homogenous and sta- 
tionnary solution ^/iq = Y^noe"*^^ *, with rig the mean 
ID density and Eq = uq U{s)ds. This solution 
is stable and can also be the ground state for small 
enough no as suggested by the Bogoliubov spectrum of 
the perturbations [26] (see below): 



Assuming that the potential scales like Uq and pos- 
sesses a single length scale a, the spectrum depends then 
only on a single dimensionless parameter |23j: 

ma^ - 

with Uq = U{s)ds ~ Uoa. For some analytical re- 
sults and for the numerics later on, we choose the soft 
core interaction, with no loss of generalities: 



U{\x-y\)^Uo9{a~\x~y\), 



(3) 



with 6{-) the Heaviside function. The Fourier transform 
of this special interaction potential is 



Uk = 2Uo 



sin(fca) 



Uo 



sin(fca) 
k ■ 



It should also be noticed that the special choice of 
the potential (|3| is purely of practical interest because 
it is easy to implement in some numerical schemes and 
can be easily used for variational estimates. Other func- 
tions whose fourier transform would be negative for a 
wavenumber domain (strictly bigger than zero), would 
show similar properties. Among them are the classi- 
cal two bodies atomic potential with strong repulsion at 
short scale and a slight attraction for large scale or a 
potential Uk choosen in such a way that the Bogoliubov 
dispersion relation matches the Landau spectrum with 
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ip{x,t) 
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my)\^dy. (4) 



Finally, we emphasize that an annular geometry can 
be simplified into a ID system by considering periodic 
boundary condition '>p{x,t) — ip{x + L,t) [L is dimen- 
sionless) and by assuming that the transverse structure 
of the solid can be neglected. We then define the energy 
and number of particles densities: 



1 

2L 



no 



A 



x+l 



L 



\'ipix,t)\'^dx. 



\iPiy)\'dy]dxX5) 



(6) 



GROUND STATE 

As shown in [53] , for low A the ground state is a super- 
fluid (without positional order) . However above a critical 
value, Ac the ground state shows a periodic modulation 
of the density in space. Although in two and three space 
dimensions the transition is first order as A increases, it 
is supercritical (second order) in one space dimension [3T|. 
The periodic structure arising from the instability can be 
analytically estimated through a variational approach for 
a fixed wavelength A at least in two regimes: close to the 
transition and for large A. 

If A ^ Ac, a weak amplitude developement of wanum- 
ber k and normalized to unity reads: 



ip{x) 



(1 -I- Ae 



ikx 



A* 



—iki 



(7) 



Minimizing the energy of such solution gives: 



2{k^+A{U2k-m)/Uo) 



(8) 



and the following wave-number selection kca — 
4.078 . . .. The amplitude for this wave number kc follows 



|Ac|2 = (A-Ac) 



-4sin(fcc) 



4fcc(fc2 + AcCos2(fcc)/3 



while for k ~ kc, \Ac\^ - \A\^ cx (fc - kef. 

In the large A limit, the density exhibits strongly non- 
linear strucures since the potential energy in ^ requires 
small ip while the mass normalization (|6| forbids to be 
small every-where. Therefore the energy minimization 
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leads to a periodic structure with zones where ■0 « bal- 
ancing zones where ip ^ 1. In the A — > cx) limit, Ref. 
[28] showed that tp ^ only in a small zone x G (—5,6) 
of the whole period (0, A). The Euler-Lagrange condition 
deduced from ^ together with ^ leads to the Hemholtz 
equation in the domain {—S,S) : —^/'{x) = ^^j. Finally 
the minimization of the energy gives 5 and the wave num- 
ber A of the periodic structure. Following this approach, 
we sketch now an estimate of the ground state for finite 
A 1. We use the trial function for a single period: 



'A f1lX\ 



(9) 



in X € [—5, S\ and zero elsewhere, that satisfies exactly 
the normalization condition (|6|. Introducing this trial 
function into the energy ([5|, we obtain: f = fi -1-^2 +^3j 
with the kinetic energy 



the self interaction of a pulse with itself 
£2 



AX 



and the nontrivial interaction of a pulse with its two near 
neighbors 



A-l-(5 



^{vfdydx. 
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This energy £3 is not zero only if A < 1 + 2(5 (and natu- 
rally, we have A > 2i5). 

The energy £ may be understood as a function of 5 
and of the periodicity length A. Then, the variation of 
total energy in the {S, A) plane shows the existence of a 
global minimum and a saddle for large enough A. As 
A decreases, the saddle and the global minimum collide, 
and we obtain a pure monotonic energy landscape in the 
(i5. A) plane, leading to both A and S to infinity as min- 
imizer. On the other hand, when A — >■ cx) the global 
minimum moves to {6, X) — > (0,1). This selection mech- 
anism holds for an infinite domain where the wavelength 
A can vary continuously. In a finite domain, A can only 
take discrete values related to the number Nc of unitary 
cells, A = L/Nc- Numerical simulations suggest also the 
existence of a large energy barrier between minimizers 
with different number of cells for large A. Thus, for 
a given domain size L, the energy, as function of (<5, A), 
is now described by a discrete set of energy functions of 
6 for each available A satisfying A = L/Nc- One has now 
to minimize each energy with respect to 5. For small A 
(typically smaller than Ac) none of these functions have 
minima. For large A on the other hand, there is a finite 
band of A for which the energy admitts a minimum as 
(5 varies. The minimization of this energy respect to 5 



provides relations among 5 and A with the wavelength A 
as a fixed parameter. To determine the global minimum 
and to avoid further algebraic difficulties, we introduce 
the new variable: z = 7r(A — l)/(5 where < z < 27r 
for our problem (in particular, z > 2-k means that the 
peaks do not interact one with another) . Minimizing the 
energy gives a relation for A = A(A, z): 



A: 



An^z 



A(A - l)2((27r - z)(cos(z) + 2) + 3sin(z)) ' 



(10) 



Fig. [T] shows 5 versus A for different values of A. The 
analytical curves are shown together with the results of 
direct numerical simulations described below. 




FIG. 1: Plot of logioS as a function of A for different A. 



The curves correspond to formula ( 10 1 while the points come 



from numerical simulations. The total size of the system is 64 
units and the range of the interaction is a = 16. The system 
displays a number of cells varying from 44 to 48 (identified 
at the right hand side of the figure). The respectives A vary 
thus from 4/3 to 16/11. 



Finally, an exponential "boundary layer" correction 
developps near x — ±(5 where the nonlinear term in eq. 
Q cannot be neglected, as noticed by [55]. In the limit 
of large A, where similarly A — > 1, the nonlinear term of 
eq. Q gives near a; = i5: 



px+l rx+l 

lim / il'(v)'^dx — Cte+lim / 



il}{yYdx « Cte+ 



The ground state is thus modified into 4'{x) + ip{x), 
where tp{x) is the trial function ([9|, and </? satisfies a 
linear Schrodinger equation: 

-lv"{x) + ^V{xMx)^0, 



where the first nontrivial term for the potential reads 
V{x) = jf^ix — S) . The solution may be com- 
puted directly in term of Bessel functions : (p{x) = 



4 



5f/^ ) . where the constant i^to one term g 



resuhs from the matching of the exponentially small 
boundary layer and the trial function. One can also ex- 
pand this solution via a WKB approxomation Lp{x) — 
^^g-x/AS(:r) ^i^-j^ 5.(2.) ^Sq{x) + ^Si{x) + .... We ob- 
tain then Sq{x) = J y/V{x)dx and Si{x) = ^log{SQ{x)) 
and therefore 



J po{x){dor:K,r:)'^ dx where 



with So(x) 



(p{x) 



"57372- 



(x-S) 
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NON-CLASSICAL MOMENT OF INERTIA IN 
SUPERFLUIDS AND SUPERSOLIDS. 

The precise estimation of the ground state is in fact 
crucial to describe the supersolid features of the model. 
Indeed, we have obtained in Refs. [241 ES], using the 
homogenization technique [29], an expression for the ef- 
fective or superfluid density matrix gf^ deduced from the 
density profile of the crystal. 

We shall in fact explore the low excited states around 
the ground state, described by the knowledge of the crys- 
tal density po{x). The change of energy for phase varia- 
tions gives: 



AE = 



1 



Po{x) 



dx 



dx. 



(11) 



and 4> is determined by minimizing A_E, that correspo, 
the Euler-Lagrange condition: 



d_ 
dx 



0. 



(12) 



under the appropriate boundary conditions. As shown 
by Leggett [3] , for a periodic po [x) under rotation, AE is 
lower than that of a rigid solid rotation, which indicates 
that superfluidity is present. 

In Refs. we have obtained an expression for 

the energy variation in three space dimensions AE = 
hn I 9ik^i4'9k4'd^Xj where g^^ is the effective or super- 
fluid density matrix. It can be explicitely expressed us- 
ing a solution of a partial differential equation in the unit 
cell V of the solid, following: 



Kx{x) is a periodic function in the interval [0, A] solution 
of d^po + {podxKx) = 0. Thus d^K^ix) = -1 + 
where c is an integration constant. The periodic bound- 
ary condition Kx{0) = Kx{X) gives c = , ^ , . Fi- 



i f 

nally, we find that in one dimension, the superfluid den- 
sity writes: 



Poix) 



dx. 



Thus, the theory of homogenization provides us an ex- 
act result for the special case of one space dimension, 
and the effective density (scalar in ID) is then a kind of 
"harmonic" average of the density [29] . 

From this formula, the non classical rotational inertia 
fraction (NCRIF) p^^/p corresponds exactly to the up- 
perbound quotient Qq proposed by Leggett [3], who also 
established the equivalence for ID systems more recently 
[30]. Therefore, the NCRIF at low speed (NCRIFq) 
reads: 



A 1 

po(x) 



dx 



(14) 



Remarks. 

1. The Schwartz inegalitY[32| and Pq{x) > gives < 
Qo < 1. 

2. For finite energy, if the ground state vanishes at 
some point, the non-classical rotational inertia does as 
well. Indeed, if at some point x^ we have po{x) ~ jcc— a;* |" 
with a > 0, then 



1 



Po{x) 



-dx ~ finite term ^ 



"dx 



and 



Qo 



1 



finite term - 



Therefore if < a < 1, Qo remains finite with e — >■ 0. 
However, as we have seen above, such a ground state 
would require an infinite amount of energy. 

RESULTS 
NCRIF in the weakly nonlinear limit. 



9ik = T^^tk - Qik 

Q^k = ^ J^Poir)VK,-VKkdr. (13) 

The vector Ki is a periodic function in the unit cell V 
that is solution of V^po + ^ ■ (Po^Ki) — 0. 

In one space dimension, we can in fact deduce the den- 



sity g^'' exactly. Indeed the formula ( 13 ) simplifies then 



The NCRIF in the limit of weak modulation may be 
computed directly from the trial function ([7|: 



2\3/2 



(1-4|A|^) 
(1 + 2|AP) 



(15) 



where |Ap is evaluated at fc = k^. As |^| 1/2, the quo- 
tient Qq vanishes, because the wavefunction ([t]) vanishes 
at some point. 
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NCRIF in the limit A ^> oo. 

For large A, since the ground state po{x) decays ex- 
ponentially, the contribution to the NCRIF ( 14 ) mainly 
comes from the large contribution of 1/pq{x) in a; e 
[(5, A/2]. That is, after using the WKB approximation: 

Qo^'-K^.e-^^^'^'-'^''\ (16) 



1-3 
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FIG. 2: (a) NCRIF as a function of v for different values of 
A. (b) NCRIFo as a function of A, the line is the curve from 
the weakly non linear analysis, see formula ( 15 1 which gives 
a good approximation up to A ~ 40. 



We will now be using numerical simulations to deduce 
the NCRI and compare it with theories by two different 
methods. First, the ground state is determined. Then, 
one can compute directly the NCRI by imposing a rota- 
tion to this ground state. On the other hand, the value 
of Qq can be calculated from the ground state solution 
po{x). 

Numerical results are obtained by minimizing the en- 
ergy ([5]) under the number of particles condition Q . We 
use therefore the Ginzburg-Landau version of the dynam- 
ics which can be interprated as the integration of the CP 
equation for imaginary time t = —it: 



dtp 



A*^+^|^-|vX^,i) / lV'(y)Pdy, (17) 



/i is the Lagrangian multiplier introduced to satisfy the 
number of particles condition. 



Imposing a rotational frequency w in a ID annular sys- 
tem amounts to consider a drift of the system at con- 
stant velocity v = uL with periodic boundary condi- 
tions. The ground state of such a system is obtained 
by minimizing: !F = £ -\- vV + p{N — no), where V — 
~Tl Ioi''P*i^)'^x^{x) — 'ip{x)dxtp*{x))dx. Consequently, 
a direct computation of the NCRIF can be performed 
numerically: 



NCRIF{v) = 1 - 



\V'{v)\ 



!,\tP{xWdx 



Fig. [2](a) shows the function NCRIF{v) for different A 
obtained by numerical minimization of J- . As expected, 
the NCRIF decreases as v increases. For large value of the 
parameter the NCRIF first become negatives and then 
show large fluctuations, indicating that complex struc- 
tures are present, such as 27r phase jumps for instance 
(similar to vortices in higher dimensions). Moreover, nu- 
merical instabilities are also enhanced by the rotation 
so that only moderate A (up to 150) values could be 
achieved with full confidence. 
The low speed limit: 



NCRIFo 



hm NCRIF{v) 



is then shown on Fig. WCb) as function of A and compared 



with the analytical quotient ( 15 ) of the weak amplitude 



modulations, showing an excellent agreement. 

On the other hand, as explained above, the NCRIFq 
can be calculated directly from the numerical solution 



Po{x), by computing the Leggett quotient Qq (14). Since 
the ground state solution is numerically more stable to 
obtain than the minimization of the rotating system, we 
are able to compute a satisfactory good estimates for Qo 
up to A of the order of 800, as shown on Fig. |3ja). 
Remarkably Qo does not depend on the wavelength of 
the periodic structure A, since all the numerical data for 
different A gather on a single curve. This is a consequence 
that the main contribution to the quotient Qo comes from 
the wide region with small values of po{x). On the other 
hand, only poor agreement is found with the asymptotic 



behavior ( 16 ). 

In Fig. [3 b) we compare this quotient Qq with the 
NCRIFq obtained by direct numerical simulation of the 
rotating system for the accessible moderate A values. It 
shows a particularly good numerical agreement between 
the two methods, as expected by the theory. 

In conclusion, we have exhibited NCRI in a ID model 
of supersolid in the context of annular geometry, us- 
ing both direct numerical simulation and analytical esti- 
mates. In particular, we have shown that the so-called 
Leggett quotient was there in full agreement with NCRI. 
C.J. acknowledges the financial support of the DGA for 
this research and this research was supported in part 
by the National Science Foundation under Grant No. 
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(a) 



800 A 




ib) 

FIG. 3: (a) The quotient Qo as a function of A, using a 
direct numerical integration of ground states po{x) obtained 
for Fig. [T| The hnes are the functions obtained from the 
theory (Tm for three different wavelength (similar notations 
as Fig. TTwith — 0.01 as a fixed parameter. Notice the 
exponential behavior in qualitative agreement with (161 . (b) 
Comparison between the numerical calculations NCRIFq and 
the Leggett's quotient Qo represented with black dots. 
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